#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "raylib.h"
#include "pub_lib.h"
#include "FreeImage.h"
#include "image.h"
#include "dpi_algo.h"

/*
 * 有了隶属度函数和规则，后面的做法就是套用之前的步骤了：
 * 一共 3 步：
 * 其实可以合并为3步：
 * 模糊处理 - 根据规则算出隶属度 - 去模糊
 *
1. 模糊处理：找对应的隶属度函, membership function
	这个函数的功能：把一个输入，映射到 [0, 1]的隶属度上
		0：隶属度最低，完全不是
		1：隶属度最高，完全就是
	对于图像处理，输入一般是灰度值 intensity,[0, 255]，
	输出是一个隶属度 float degree, [0, 1]
	typedef float (membership*)(BYTE val);

2.  根据事物逻辑和推断方法，计算出最终的隶属度 Q
	这个Q 只和对应的期望值 V 相关

3. 去模糊：得到一个“干脆的”标量输出（比如重心公式)
	out = $\frac{\sum_{v=1}^K  V_i * Q(V_i)} {\sum_{v=1}^K Q(V_i)}$

	Vi: 期望值，一共有K 种
	Q(Vi): 对于当前期望值Vi，对应的隶属度，值在 [0, 1] 范围内
	分子 Numerator  ：期望值 * 对应的隶属度之和，一共 K 对乘积
	分母 Denominator：K种隶属度之和
	重心坐标 = Numerator / Denominator
 * */

/* ch3.8.4 使用模糊集合进行灰度变换
 * 输入 z ：当前像素的灰度值
 *	最终计算结果只和当前像素的灰度值相关
 *
 * 3个期望值：Vd = 0, Vg = 127, Vb = 255
 *
 * 1. 确定隶属度函数: u_dark, u_gray, u_bright
 *	插图 3.55，图c
 *
 * 2. 根据事物逻辑和推断方法，计算出最终的隶属度 Q
 *	如果一个像素是暗的，则它的隶属度 Q(Vd) = dark(z)
 *	如果一个像素是灰的，则它的隶属度 Q(Vg) = gray(z)
 *	如果一个像素是亮的，则它的隶属度 Q(Vb) = bright(z)
 *
 * 3. 去模糊：得到一个“干脆的”标量输出（比如重心公式)
 * 3个期望值 Vd, Vg, Vb，所以 K = 3
 *	重心公式简化为：
 *		Q(v): dark(z), gray(z), bright(z)
 *		v:		Vd,		Vg,		Vb
 *		分子 = dark(z) * Vd + gray(z) * Vg + bright(z) * Vb;
		分母 = dark(z) + gray(z) + bright(z);
		output = sum / div;
 * */

/* 输入一个灰度值 val
 * 输出归一化的隶属度 [0, 1]
 * */
typedef float (*u_func)(BYTE val);

/* Trapezoid，梯形
 * a = c = 0
 * b = 63, b + d = 127, d = 64
 * */
static float u_dark(BYTE val)
{
	float ret = 0.0f;
	if(val < 63)
		ret = 1.0f;
	else if(val >= 63 && val < 127)
		ret = 1.0f - (val - 63) / 64.0f;
	else
		val = 0.0f;
	return ret;
}

/* Triangle，三角形
 * a = 127
 * a - b = 63, b = 64
 * a + c = 192, c = 65
 * */
static float u_gray(BYTE val)
{
	float ret = 0.0f;
	if(val < 63)
		ret = 0.0f;
	else if(val >= 63 && val < 127)
		ret = 1.0f - (127 - val) / 64.0f;
	else if(val >= 127 && val <= 192)
		ret = 1.0f - (val - 127) / 65.0f;
	else
		ret = 0.0f;
	return ret;
}

/* Sigma
 * a = 192
 * a - b = 127, b = 65
 * */
static float u_bright(BYTE val)
{
	float ret = 0.0f;
	if(val < 127)
		ret = 0.0f;
	else if(val >= 127 && val <= 192)
		ret = 1.0f - (192 - val) / 65.0f;
	else
		ret = 1.0f;
	return ret;
}

//cal center of gravity of z
BYTE cal_center_gravity(u_func dark, u_func gray, u_func bright,\
		BYTE Vd, BYTE Vg, BYTE Vb, BYTE z)
{
	float sum, div, ret;
	sum = dark(z) * Vd + gray(z) * Vg + bright(z) * Vb;
	div = dark(z) + gray(z) + bright(z);
	ret = sum / div;
	return (BYTE)ret;
}

/* figure 3.54 */
void fuzzy_intensity_hence(FIBITMAP* ori, FIBITMAP* chg)
{
	BYTE Vd = 0;
	BYTE Vg = 127;
	BYTE Vb = 255;

	unsigned int x, y;
	BYTE value;

	unsigned int w, h;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);

	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(ori, x, y, &value);
			value = cal_center_gravity(u_dark, u_gray, u_bright,\
				Vd, Vg, Vb, value);
			FreeImage_SetPixelIndex(chg, x, y, &value);
		}
	}
}


/* ch3.8.5 使用模糊集合进行空间滤波
 * 输入：当前像素的灰度值 z
 * 3套规则：
 *	公共规则1：ZE，当前像素是 "0" 的隶属度，图3.57a，高斯分布
 *
 *	计算 ”黑色“隶属度
 *	规则1：BL，当前像素是 "黑色" 的隶属度，图3.57b，梯形
 *	规则2：图3.58所示，根据 3*3 邻域，得到中心点的 WH 白色 隶属度
 *		BL = 1 - WH
 *
 *	计算 “白色”的隶属度：
 *	规则1：WH，当前像素是”白色“的隶属度，图3.57b，Sigma形
 *
 *	规则2：图3.58所示，根据 3*3 邻域，得到中心点的 WH 白色隶属度
 *		本质：中心点和周围点的灰度差接近0，则中心点为白色，很亮
 *		中心点和周围点灰度差越大，则中心点是黑色的隶属度越大
 *		这样可以确定边界，轮廓
 *
 *	256 个可能的期望值，[0, 255]，所以 K = 256
 *
 * 1. 隶属度函数: ZE, BL, WH
 *	插图 3.57
 *
 * 2. 根据事物逻辑和推断方法，计算出最终的隶属度 Q
 *
 * 先应用规则2， 图3.58所示，根据 3*3 邻域，得到中心点的 WH 白色隶属度
 *
 * if d2 and d6 is 0, then z5 is white
 *	这里 0, white 都是模糊的
 *
 *  "0" 对应的隶属度函数: ZE
 *  对于这个邻域 d2, d6：因为需要同时满足, 所以去交集，min
 *  当前点
 *		d2 = z2 - center
 *		v2 = ZE(d2) = fuzzy_set_zero(d2)；根据灰度差计算属于0的隶属度
 *		同理
 *		d6 = z6 - center
 *		v6 = ZE(d6)
 *
 *		WH1 = min(v2, v6) d2 and d6 is 0, then white
 *
 * 图3.58 的规则还包括以下的：
 *    if d2 and d6 is 0, then z5 is white, wh1,
 *	  if d6 and d8 is 0, then z5 is white, wh2,
 *	  if d8 and d4 is 0, then z5 is white, wh3
 *	  if d4 and d2 is 0, then z5 is white, wh4
 *		计算出4个邻域的隶属度后，满足这4个条件的任何一个，中心点都是白色
 *		当前点并集，最大值，max
 *		wh = max(wh1, wh2, wh3, wh4)
 *
 *	应用规则2，最终计算出当前点的”白色“的隶属度： wh
 *		当前点是”黑色“的隶属度：bl = 1 - wh
 *	对应的实现：fuzzy_calc_wh
 *
 * 3. 去模糊：得到一个“干脆的”标量输出（比如重心公式)
 *
 *	黑色的隶属度是 BL
 *	白色的隶属度是 WH
 *
 * k: 256种可能的值，0 - 255, K = 256
 * Vi: 期望值，即 0 - 255
 * Q(Vi): 对于每一个期望值，最终的隶属度
 *
 * 去模糊计算：fuzzy_calc_gravity
 *
 * 规则1：BL，当前像素是 "黑色" 的隶属度，图3.57b，梯形
 *	根据规则2的计算，得到当前像素，”黑色“隶属度是bl
 *  还需要满足规则1，fuzzy_set_bl() 计算出来的隶属度
 *  需要同时满足2个条件，所以用交集，min
 *		black = min(规则1，bl)
 *
 *	注意：规则2是根据图片”邻域“计算得到的隶属度
 *		规则1是根据”灰度值“得到的隶属度
 *		需要同时满足这两个规则，才能说这个像素是”黑色“的
 *
 *	同理：
 *		white = min(规则1，wh)
 *
 *  一个像素，模糊处理后，
 *		”黑色“的隶属度 black
 *		”白色“的隶属度 white
 * 最终的隶属度，取最大值，并集，max
 * 对于同一个输入z（中心点的灰度值），只能有一个隶属度，取最大的
 *		Q = max(black, white)
 *
 * 伪代码
 *
 *  bl = 1 - wh
 *	for(i = 0; i < 256; i++) {
 *		black = min(规则1，bl)
 *		white = min(规则1，wh)
 *		Q = max(black, white)
 *		sum += Q * i
 *		div += Q
 *	}
 *
 * 	最终灰度 = 重心值 = sum / div
 *
 * 	相比第一个例子：
 *	第1个例子： 只有1个规则，所以无需取并集或交集
 *		隶属度只和当前像素的灰度值相关
 *		期望值是常数，且只有3个期望值，容易计算
 *
 *	第2个例子：3个规则，并且第2个规则比较繁琐
 *		隶属度不仅和当前像素相关，还和邻域相关
 *		期望值也不是常数，也是模糊的
 *		期望值被离散成 [0, 255] 256个数字
 *		所以计算复杂得多
 * */

/*
 * ZE: Zero
 * ZE(z) = exp(K * z * z)
 * 从插图 3.57(b) 看出，灰度值50时，隶属度接近0，这里取值 0.1
 * Z = 50, Ze(50) = 0.1, so K = -0.000921
 */ 
float fuzzy_set_zero(int z)
{
	double K = -0.000921;
	double ret;
	ret = exp(K * z * z);
	return (float)ret;
}

/*
 * BL: Black
 * 200这个数字从何而来：
 * 从插图 3.57(b) 看出: 灰度值200时，隶属度为0
 * 灰度值0 时，隶属度为 1
*/
float fuzzy_set_bl(BYTE z)
{
	float ret = 0.0f;
	if(z < 200)
		ret = (200 - z) / 200.0f;
	return ret;
}

//WH: White
float fuzzy_set_wh(BYTE z)
{
	float ret = 0.0f;
	if(z > 50)
		ret = (z - 50) / 200.0f;
	return ret;
}

typedef struct {
	float BL[256];
	float WH[256];
	BYTE area[9];
}fuzzy_extra;

float fuzzy_calc_wh(BYTE* area)
{
	int d2, d4, d6, d8;
	float v2, v4, v6, v8;
	float wh1, wh2, wh3, wh4;
	int center;

	center = area[4];
	//figure 2.56 (b)
	d2 = area[1] - center;
	d4 = area[3] - center;
	d6 = area[5] - center;
	d8 = area[7] = center;

	//figure 3.57 (a)，计算 dx 的隶属度
	v2 = fuzzy_set_zero(d2);
	v4 = fuzzy_set_zero(d4);
	v6 = fuzzy_set_zero(d6);
	v8 = fuzzy_set_zero(d8);

/*    if d2 and d6 is 0, then z5 is white
 *	  if d6 and d8 is 0, then z5 is white
 *	  if d8 and d4 is 0, then z5 is white
 *	  if d4 and d2 is 0, then z5 is white
 *		因为每一个条件中是and，所以取交集，用 min
 * */
	wh1 = min(v2, v6);
	wh2 = min(v6, v8);
	wh3 = min(v8, v4);
	wh4 = min(v4, v2);
/*
 * 这4个条件，满足任何一个，中心点都是白色，所以取并集，用 max
* */
	float wh5, wh6, wh;
	wh5 = max(wh1, wh2);
	wh6 = max(wh3, wh4);
	wh = max(wh5, wh6);
	return wh;
}

BYTE fuzzy_calc_gravity(float wh, float* pbl, float* pwh)
{
	float bl, sum, div, Q;
	float black, white;
	int i;
	bl = 1.0f - wh;

	sum = 0.0f;
	div = 0.0f;
	for(i = 0; i < 256; i++) {
		//1个输入，必须同时满足这2个规则，所以取交集, min
		black = min(pbl[i], bl);
		white = min(pwh[i], wh);
		//哪个隶属度值大，取哪个，所以用并集, max
		Q = max(black, white);
		sum += i * Q;
		div += Q;
	}
	
	float ret;
	ret = sum / div;
	return (BYTE)ret;
}

static void fuzzy_filter_example(unsigned int index, int* sum, unsigned char k,\
	int dim, BYTE val, BYTE* new_val, void* extra)
{
	BYTE *area;
	fuzzy_extra* ptr;
	int end = dim * dim;

	//fill area[9]
	ptr = (fuzzy_extra*)extra;
	area = ptr->area;
	area[k - 1] = val;

	float wh;
	BYTE gravity;
	if (k == end) {
		wh = fuzzy_calc_wh(area);
		gravity = fuzzy_calc_gravity(wh, ptr->BL, ptr->WH);
		*new_val = gravity;
	}
}

/* 3.8.5 */
void fuzzy_spatial_filter_hence(FIBITMAP* ori, FIBITMAP* chg)
{
	int i;
	fuzzy_extra extra;

	//init BL, WH
	for(i = 0; i < 256; i++) {
		extra.BL[i] = fuzzy_set_bl(i);
		extra.WH[i] = fuzzy_set_wh(i);
	}

	/* 3 * 3 area */
	imgage_filter_common(ori, chg, 3, fuzzy_filter_example, (void*)&extra);
}

